function elements = calculateLimitStage(elements,model)
ne = model.countElements();
for i=1:ne
    e = model.getElement(i-1);
    face = e.getFace();
    ue = face.getUHatLocal();
    xi = e.getGaussPoint(0);
    B = zeros(3,6);
    B = e.computeB(face,B,xi);
    %check if the crack was initialized
    if elements(i).isCrack==0 && elements(i).isActivated==1
        sigma = elements(i).stress;
        %perform spectral decomposition
        oxx = sigma(1);
        oyy = sigma(2);
        oxy = sigma(3);
        t1 = (oxx + oyy) / 2;
        t2 = sqrt(((oxx - oyy)/2)^2 + oxy^2);
        o1 = t1 + t2;
        if o1>elements(i).ft
            disp(['element ' num2str(i) ' crack initiated']);
            elements(i).isCrack = 1;
            %check the solitary node
            if abs(oxy) < 1e-12
                if (oxx > oyy)
                    n1 = 1;
                    n2 = 0;
                else
                    n1 = 0;
                    n2 = 1;
                end
            else
                g = sqrt((oxx - o1)^2 + oxy^2);
                n1 = -oxy / g;
                n2 = (oxx - o1) / g;
            end
            N = face.getNodes();
            cen = face.computeCenter();
            cx = cen(1);
            cy = cen(2);
            d = -n1*cx - n2*cy;
            for j = 1:3
                X = N(j).getPosition();
                f(j) = n1*X(1) + n2*X(2) + d;
            end
            if f(1) * f(2) > 0
                PosDomain = [3];
                NegDomain = [1 2];
            else
                if f(1) * f(3) > 0
                    PosDomain = [2];
                    NegDomain = [1 3];
                else
                    PosDomain = [1];
                    NegDomain = [2 3];
                end
            end
            %compute the crack normal direction
            posNode = N(PosDomain(1)).getPosition();
            delta(1) = posNode(1) - cx;
            delta(2) = posNode(2) - cy;
            if(n1*delta(1)+n2*delta(2) < 0)
                n1 = -n1;
                n2 = -n2;
            end
            elements(i).crack_normal = [n1 n2]';
            %compute intersection for post-processing
            xi1 = posNode(1);
            yi1 = posNode(2);
            negNode1 = N(NegDomain(1)).getPosition();
            negNode2 = N(NegDomain(2)).getPosition();
            xi2 = negNode1(1);
            yi2 = negNode1(2);
            xi3 = negNode2(1);
            yi3 = negNode2(2);
            a1 = yi2 - yi1;
            b1 = xi1 - xi2;
            c1 = yi1 * (xi2 - xi1) - xi1 * (yi2 - yi1);
            a2 = yi3 - yi1;
            b2 = xi1 - xi3;
            c2 = yi1 * (xi3 - xi1) - xi1 * (yi3 - yi1);
            d = -n1*cx - n2*cy;
            X1 = computeIntersection(n1, n2, d, a1, b1, c1);
            X2 = computeIntersection(n1, n2, d, a2, b2, c2);
            elements(i).X1 = X1;
            elements(i).X2 = X2;
            %compute H
            H = zeros(6,2);
            H(2*PosDomain-1,:) = [n1 -n2];
            H(2*PosDomain,:) = [n2 n1];
            elements(i).H = H;
            %compute P
            P = [n1^2 -n1*n2;
                n2^2 n1*n2;
                2*n1*n2 n1^2-n2^2]; %TODO: check the accuracy of P
            C = elements(i).constitutiveMatrix;
            %compute A
            A = P'*C*B;
            elements(i).A = A;
%             %postulate Dnn to ensure traction consistency and unique solution
%             t0 = A*ue;
%             ft = elements(i).ft;
%             Dnn = ft^2/(t0'*(Dhat\t0));
%             elements(i).Dnn = Dnn;
%             elements(i).Dss = Dnn;
%             Dhat = [Dnn 0;0 Dss];
%             %g0 = sqrt(Dnn*(t0'*(Dhat\t0)));
%             %elements(i).ttilde = g0;
%             %compute v0
%             v0 = Dnn*(Dhat\t0)/g0;
%             %validate v0
%             assert(norm(t0'*v0-g0)<1e-9,['element' num2str(i) ' violate condition t0''*v0-g0=0']);
%             assert(norm(v0'*Dhat*v0-Dnn)<1e-9,['element' num2str(i) ' violate condition v0''*Dhat*v0-Dnn=0']);
%             %compute f0
%             f0 = Dhat*v0/g0;
%             elements(i).f = f0;
            %set initial state for element
            elements(i).isOpening = 1;
            elements(i).isSticking = 0;
            elements(i).isSliding = 0;
            elements(i).isLoading = 1;
            %compute initial crack opening & gamma
            [e,gamma] = computeInitialOpening(elements(i),ue,[1e-4 1e-4]');
%             [e,gamma] = computeGammaOpening(elements(i),ue,[1e-4 1e-4]')
            %validate e and gamma
            Dnn = elements(i).Dnn;
            Dss = elements(i).Dss;
            Dhat = [Dnn 0;0 Dss];
            etilde = computeEquivalentSeparation(elements(i),e);
            assert(abs(gamma-computeDamageFunction(elements(i),etilde))<1e-14,...
                ['element ' num2str(i) ': error gamma != F(etilde(e))']);
%             norm((Dhat+gamma*A*H)*e-gamma*A*ue)
            assert(norm((Dhat+gamma*A*H)*e-gamma*A*ue)<1e-14,...
                    ['element ' num2str(i) ' violate traction equilibrium']);
            elements(i).e = e;
            elements(i).gamma = gamma;
            %update t
            t = Dhat*e/gamma;
            elements(i).t = t;
            %update f
            dFdetilde = computeDamageFunctionDerivative(elements(i),etilde);
            elements(i).f = dFdetilde/etilde*gamma/Dnn*t;
        end
    end
end
end

% function [e,gamma] = computeInitialOpening(element,d,e0)
% e = e0;
% etilde = computeEquivalentSeparation(element,e);
% gamma = computeDamageFunction(element,etilde);
% A = element.A;
% H = element.H;
% Dnn = element.Dnn;
% Dss = element.Dss;
% Dhat = [Dnn 0;0 Dss];
% rhs = gamma*A*d - (Dhat + gamma*A*H)*e;
% f0 = element.f;
% conv_criteria = 0;
% % element.ef
% while conv_criteria == 0
%     lhs = Dhat + gamma*A*H + A*(H*e-d)*f0';
%     de = lhs\rhs;
%     e = e+de;
%     etilde = computeEquivalentSeparation(element,e);
%     gamma = computeDamageFunction(element,etilde);
%     rhs = gamma*A*d - (Dhat + gamma*A*H)*e;
%     if norm(rhs) < 1e-14
%         conv_criteria = 1;
%     end
% %     gamma*(A*d-A*H*e)-Dhat*e
% %     rhs
% end
% end


% function [e,gamma] = computeInitialOpening(element,d,e0)
% e = e0;
% etilde = computeEquivalentSeparation(element,e);
% gammas = 1/computeDamageFunction(element,etilde);
% A = element.A;
% H = element.H;
% Dnn = element.Dnn;
% Dss = element.Dss;
% Dhat = [Dnn 0;0 Dss];
% rhs = A*d - (gammas*Dhat + A*H)*e;
% conv_criteria = 0;
% while conv_criteria == 0
%     dFdetilde = computeDamageFunctionDerivative(element,etilde);
%     f = dFdetilde*Dhat*e/Dnn/etilde;
%     F = computeDamageFunction(element,etilde);
%     lhs = gammas*Dhat + A*H - Dhat*e*f'/F^2;
%     de = lhs\rhs;
%     e = e+de;
%     etilde = computeEquivalentSeparation(element,e);
%     gammas = 1/computeDamageFunction(element,etilde);
%     rhs = A*d - (gammas*Dhat + A*H)*e;
%     if norm(rhs) < 1e-14
%         conv_criteria = 1;
%     end
% end
% gamma = 1/gammas;
% end

function [e,gamma] = computeInitialOpening(element,d,e0)
e = e0;
etilde = computeEquivalentSeparation(element,e);
gammas = 1/computeDamageFunction(element,etilde);
A = element.A;
H = element.H;
Dnn = element.Dnn;
Dss = element.Dss;
Dhat = [Dnn 0;0 Dss];
conv_criteria = 0;
while conv_criteria == 0
    e = (gammas*Dhat + A*H)\(A*d);
    etilde = computeEquivalentSeparation(element,e);
    gammas = 1/computeDamageFunction(element,etilde);
    rhs = A*d - (gammas*Dhat + A*H)*e;
    if norm(rhs) < 1e-14
        conv_criteria = 1;
    end
end
gamma = 1/gammas;
end

% function [e,gamma] = computeInitialOpening(element,d,e0)
% A = element.A;
% H = element.H;
% Dnn = element.Dnn;
% Dss = element.Dss;
% Dhat = [Dnn 0;0 Dss];
% e = e0;
% etilde = computeEquivalentSeparation(element,e);
% gamma = computeDamageFunction(element,etilde);
% to = A*(d - H*e);
% te = Dhat*e/gamma;
% rhs = gamma*(to - te);
% conv_criteria = 0;
% 
% while conv_criteria == 0
%     dFdetilde = computeDamageFunctionDerivative(element,etilde);
%     f = dFdetilde*Dhat*e/Dnn/etilde;
%     lhs = Dhat + gamma*A*H - to*f';
%     de = lhs\rhs;
%     e = e+de;
%     etilde = computeEquivalentSeparation(element,e);
%     gamma = computeDamageFunction(element,etilde);
%     to = A*(d - H*e);
%     te = Dhat*e/gamma;
%     rhs = gamma*(to - te);
%     if norm(rhs) < 1e-14
%         conv_criteria = 1;
%     end
% end
% end
